library(gridExtra)
source('loader.R')

getfilename <- function(filename)
	paste('../tex/', filename, '.tex', sep='')

output <- function(filename, width)
	pdf(paste('../tex/', filename, '.pdf', sep=''),
		paper='special', width=8, height=8/width)


## modify ggplot2 graphics to make more publishing-like
theme_set(
	theme_gray(base_size=12, base_family="") %+replace% 
        theme(
        	axis.text = element_text(colour="black"),
            panel.grid.minor = element_blank(),
            panel.grid.major = element_blank(),
            panel.background = element_rect(fill="white"),
            legend.key = element_blank()
    )
)


X <- read('paiva')


output('paiva-data', 1.1)
p1 <- showdata(X, -1, 15)
p2 <- showdata(X, -1, 0)
grid.arrange(p1, p2, nrow=2)
dev.off()


states <- c('0')
for(s in states) {
	output(paste('paiva-lambdas-hist-', s, sep=''), 2)
	lambdas_hist(Xstate(X, s), c(-0.2,1), c(-4,+8))
	dev.off()
}






states <- c('0','00','000','0000')
output('paiva-lambdas-norm', 4)
lambdas_fit('paiva', states, pnorm, c(-0.1,0.4), c(0,0.0028))
dev.off()
output('paiva-lambdas-bhp', 4)
lambdas_fit('paiva', states, pbhp, c(0.2,0.7), c(0,0.0028))
dev.off()






p1 <- best_lambdas('paiva', '0', 12, test.distance, min)
p1$layers <- c(geom_hline(yintercept=0.00075,colour="gray",linetype='dashed'), p1$layers)

p2. <- best_lambdas('paiva', '1', 12, test.distance, min) + theme(legend.position="none")
p2bottom <- p2. + scale_y_continuous(breaks=c(0.00025,0.00075)) + coord_cartesian(xlim=c(-0.05,0.3), ylim=c(0,0.001)) + theme(plot.title=element_blank(),axis.title.y=element_blank())
p2bottom$layers <- c(geom_hline(yintercept=0.00075,colour="gray",linetype='dashed'), p2bottom$layers)

p2top <- p2. + theme(axis.title.x=element_blank())
p2top$layers <- c(geom_rect(size=0.3, fill='#e1e1e1', color='#e1e1e1', xmin=-0.05, xmax=0.3, ymin=0, ymax=0.001), p2top$layers)

library(gtable)  # top and bottom
g <- gtable:::rbind_gtable(ggplotGrob(p2top), ggplotGrob(p2bottom), "first")
panels <- g$layout$t[grep("panel", g$layout$name)]
g$heights[panels] <- lapply(c(4,1), unit, "null")  # different heights
p2 <- gtable_add_rows(g, unit(-1,"cm"), pos=nrow(ggplotGrob(p2top)))
p2 <- arrangeGrob(p2)


# truque para legenda comum
g_legend<-function(p){
	tmp <- ggplotGrob(p)
	leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
	dev.off()
	legend <- tmp$grobs[[leg]]
	return(legend)
}
legend <- g_legend(p1)
lwidth <- sum(legend$width[2])

output('paiva-mem-fit', 2.30)
grid.arrange(
	p1+theme(legend.position="none"),
	p2,
	legend,
	widths=unit.c(unit(0.56,'npc')-lwidth, unit(0.56,'npc')-lwidth, lwidth),
	nrow=1)
dev.off()



print("Fisher tests !")
print("state1 state2 pvalue OR")

for(m in 1:5) {
	t <- memmarkov(X, m, FALSE)
	n <- 2^(m-1)
	for(i in 1:n) {
		j <- i+n
		f <- fisher.test(matrix(c(t[,i], t[,j]), nrow=2))
		pvalue <- f$p.value
		or <- f$estimate
		si <- getsymbol(i-1, m)
		sj <- getsymbol(j-1, m)
		print(paste(format(si,width=6), format(sj,width=6), format(round(pvalue,2),width=6), round(or,2)))
	}
}


# global warming - interesting:

#states <- c('0','00','0000','000000')
states <- c('0','000000')
for(s in states) {
	output(paste('paiva-lambdas-time-', s, sep=''), 2.5)
	print(lambdas_time('paiva', s))
	dev.off()
}

